extern "C" __device__ 
double signed_triangle_area(int ax, int ay, int bx, int by, int cx, int cy) 
{
    return 0.5 * ((by-ay)*(bx+ax) + (cy-by)*(cx+bx) + (ay-cy)*(ax+cx));
}

extern "C" __device__ 
double min3(int a, int b, int c) 
{
    return max(min(min(a, b), c), 0);
}
extern "C" __device__ 
double max3(int a, int b, int c, int bound) 
{
    return min(max(max(a, b), c), bound);
}

extern "C" __global__ 
void triangle(int ax, int ay, float az, int bx, int by, float bz, int cx, int cy, float cz, int width, int height, uchar4* framebuffer, float* z_buffer, uchar4 color) {
    int bbminx = min3(ax, bx, cx); // bounding box for the triangle
    int bbminy = min3(ay, by, cy); // defined by its top left and bottom right corners
    int bbmaxx = max3(ax, bx, cx, width - 1);
    int bbmaxy = max3(ay, by, cy, height - 1);
    double total_area = signed_triangle_area(ax, ay, bx, by, cx, cy);
    if (total_area < 1) return;

    for (int x=bbminx; x<=bbmaxx; x++) {
        for (int y=bbminy; y<=bbmaxy; y++) {
            double alpha = signed_triangle_area(x, y, bx, by, cx, cy) / total_area;
            double beta  = signed_triangle_area(x, y, cx, cy, ax, ay) / total_area;
            double gamma = signed_triangle_area(x, y, ax, ay, bx, by) / total_area;
            if (alpha >= 0 && beta >= 0 && gamma >= 0) {
                float z = alpha * az + beta * bz + gamma * cz;
                int fb_idx = x + y * width;
                if (z > z_buffer[fb_idx])
                {
                    z_buffer[fb_idx] = z;
                    framebuffer[fb_idx] = color;
                }    
            }
        }
    }
}

extern "C" __device__ 
int simple_rand(int* state) {
    // XOR-shift 算法
    int x = *state;
    x ^= x << 13;
    x ^= x >> 17;
    x ^= x << 5;
    *state = x;
    return x;
}

extern "C" __device__ 
void mat4_mul_vec4(float* mat4, float* vec4) {
    float x0 = mat4[0] * vec4[0] + mat4[1] * vec4[1] + mat4[2] * vec4[2] + mat4[3] * vec4[3];
    float x1 = mat4[4] * vec4[0] + mat4[5] * vec4[1] + mat4[6] * vec4[2] + mat4[7] * vec4[3];
    float x2 = mat4[8] * vec4[0] + mat4[9] * vec4[1] + mat4[10] * vec4[2] + mat4[11] * vec4[3];
    float x3 = mat4[12] * vec4[0] + mat4[13] * vec4[1] + mat4[14] * vec4[2] + mat4[15] * vec4[3];
    vec4[0] = x0;
    vec4[1] = x1;
    vec4[2] = x2;
    vec4[3] = x3;
}

extern "C" __device__ 
void vec4_div_w(float* vec4) {
    vec4[0] = vec4[0] / vec4[3];
    vec4[1] = vec4[1] / vec4[3];
    vec4[2] = vec4[2] / vec4[3];
    vec4[3] = 1;
}

extern "C" __device__ 
void cross(float* a, float* b, float* ret) {
    ret[0] = a[1] * b[2] - a[2] * b[1];
    ret[1] = a[2] * b[0] - a[0] * b[2];
    ret[2] = a[0] * b[1] - a[1] * b[0]; 
}
extern "C" __device__ 
void norm(float* x) {
    float len = x[0] * x[0] + x[1] * x[1] + x[2] * x[2];
    len = sqrt(len);
    x[0] = x[0] / len;
    x[1] = x[1] / len;
    x[2] = x[2] / len;
}
extern "C" __device__ 
float dot(float* a, float* b) {
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
extern "C" __device__ 
void reflect(float* n, float* l, float* ret) {
    float cos2 = dot(n, l) * 2;
    ret[0] = n[0] * cos2;
    ret[1] = n[1] * cos2;
    ret[2] = n[2] * cos2;
    ret[0] -= l[0];
    ret[1] -= l[1];
    ret[2] -= l[2];
}


extern "C" __global__ 
void render_model_kernel(
    float* model_verts, 
    unsigned int* model_face_vrts,
    int face_num,
    char* frame_buffer, 
    float* z_buffer, 
    int width, 
    int height,
    float* model_view,
    float* perspective,
    float* viewport)
{
    int face_idx = blockIdx.x * blockDim.x + threadIdx.x; 
    if (face_idx < face_num) 
    {
        unsigned int v1_idx = model_face_vrts[face_idx * 3];
        unsigned int v2_idx = model_face_vrts[face_idx * 3 + 1];
        unsigned int v3_idx = model_face_vrts[face_idx * 3 + 2];

        float v1[4];
        v1[0] = model_verts[v1_idx * 4];
        v1[1] = model_verts[v1_idx * 4 + 1];
        v1[2] = model_verts[v1_idx * 4 + 2];
        v1[3] = 1;

        float v2[4];
        v2[0] = model_verts[v2_idx * 4];
        v2[1] = model_verts[v2_idx * 4 + 1];
        v2[2] = model_verts[v2_idx * 4 + 2];
        v2[3] = 1;

        float v3[4];
        v3[0] = model_verts[v3_idx * 4];
        v3[1] = model_verts[v3_idx * 4 + 1];
        v3[2] = model_verts[v3_idx * 4 + 2];
        v3[3] = 1;

        mat4_mul_vec4(model_view, v1);
        mat4_mul_vec4(model_view, v2);
        mat4_mul_vec4(model_view, v3);

        float ambient = 0.3;

        float l[3];
        l[0] = 1;
        l[1] = 1;
        l[2] = 1;

        float v1_v2[3];
        v1_v2[0] = v2[0] - v1[0];
        v1_v2[1] = v2[1] - v1[1];
        v1_v2[2] = v2[2] - v1[2];

        float v1_v3[3];
        v1_v3[0] = v3[0] - v1[0];
        v1_v3[1] = v3[1] - v1[1];
        v1_v3[2] = v3[2] - v1[2];

        float n[3];
        cross(v1_v2, v1_v3, n);
        norm(n);
        norm(l);
        float diffuse = dot(l, n);

        float r[3];
        reflect(n, l, r);
        float r_dot_v = r[2] > 0 ? r[2] : 0;
        float specular = pow(r_dot_v, 35);

        uchar4 color;
        float coe = ambient + 0.4 * diffuse + 0.9 * specular;
        coe = coe > 1 ? 1 : coe;
        color.x = 255 * coe;
        color.y = 255 * coe;
        color.z = 255 * coe;
        color.w = 255;

        mat4_mul_vec4(perspective, v1);
        mat4_mul_vec4(perspective, v2);
        mat4_mul_vec4(perspective, v3);

        vec4_div_w(v1);
        vec4_div_w(v2);
        vec4_div_w(v3);

        mat4_mul_vec4(viewport, v1);
        mat4_mul_vec4(viewport, v2);
        mat4_mul_vec4(viewport, v3);

        triangle(
            v1[0], v1[1], v1[2],
            v2[0], v2[1], v2[2],
            v3[0], v3[1], v3[2],
            width, height,
            (uchar4 *)frame_buffer,
            z_buffer,
            color
        );
    }
}